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ABSTRACT 

We perform a cosmological simulation in order to model the growth and evolution 
of Population III (Pop III) stellar systems in a range of host minihalo environments. 
A Pop III multiple system forms in each of the ten minihaloes, and the overall mass 
function is top-heavy compared to the currently observed initial mass function in the 
Milky Way. Using a sink particle to represent each growing protostar, we examine the 
binary characteristics of the multiple systems, resolving orbits on scales as small as 20 
AU. We find a binary fraction of ~ 36%, with semi-major axes as large as 3000 AU. 
The distribution of orbital periods is slightly peaked at < 900 yr, while the distribution 
of mass ratios is relatively flat. Of all sink particles formed within the ten minihaloes, 
~ 50% are lost to mergers with larger sinks, and ~ 50% of the remaining sinks are 
ejected from their star-forming disks. The large binary fraction may have important 
implications for Pop III evolution and nucleosynthesis, as well as the final fate of the 
first stars. 
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theory - first stars - early Universe 


1 INTRODUCTION 

The period during which the first stars arose marked a piv- 
otal turning point in the evolution of the Universe. The first 
stars, also known as Population III (Pop III), are believed 
to have formed at z > 20 within dark matter (DM) mini- 
haloes of 10 6 Mg (e.g. Haiman et al. 1996; Tegmark et al. 
1997; Yoshida et al. 2003). Afterwards they generated ion- 
izing photons that contributed to the reionization of 
the Universe (e.g. Kitayama et al. 2004; Sokasian et al. 
2004; Whalen et al. 2004; Alvarez et al. 2006; Johnson et al. 
2007), while those Pop III stars that ended their 
lives as supernovae (SNe) contributed to the metallic- 
ity of the intergalactic medium (IGM; Madau et al. 2001; 
Mori et al. 2002; Bromm et al. 2003; Wada & Venkatesan 
2003; Norman et al. 2004; Tornatore et al. 2007; Greif et al. 
2007, 2010; Wise & Abel 2008; Maio et al. 2011; recently re- 
viewed in Karlsson et al. 2012). This metal-enrichment and 
growing radiation background provided by the first stars de- 
termined the environment in which later stellar generations 
formed. 

The mass of a given Pop III star plays the central role in 
determining how the star will influence its surroundings. For 
instance, more massive stars are more efficient in produc- 
ing ionizing photons, per unit stellar mass, than low-mass 
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stars (Bromm et al. 2001; Schaerer 2002). Furthermore, non- 
rotating primordial stars with mass 140 M@ < M, < 260 Mg 
are believed to have exploded as pair-instability supernovae 
(PISNe; Heger & Woosley 2002), releasing the entirety of 
their metal content into the IGM and surrounding haloes, 
while stars within the range 15 Mq < M» < 40 M© ended 
their lives as core-collapse SNe. On the other hand, non- 
rotating Pop III stars with main sequence masses in the 
range 40 M© < M* < 140 M© or M* > 260 M© are expected 
to collapse directly into black holes (BHs) , thus contributing 
no metals to their surroundings. 

BH remnants of Pop III stars may emit copious amounts 
of radiation during subsequent mass accretion, however. 
While this radiation will heat and ionize its local surround- 
ings, the X-ray component has a long mean-free path that 
allows it to ionize the IGM and minihalo gas as far as 
10 — 100 kpc away. This contributed to reionization, while 
the boosted electron fraction can catalyze the formation 
of additional H 2 molecules, mildly enhancing the global 
cooling and star-formation rate (e.g., Madau et al. 2004; 
Ricotti & Ostriker 2004; Kuhlen & Madau 2005; Haiman 
2011; Jeon et al. 2012). Recent work has demonstrated that 
a BH remnant which has a stellar binary companion, giv- 
ing rise to a high-mass X-ray binary (HMXB), will exert 
stronger feedback on its surroundings than an isolated BH 
(e.g., Mirabel et al. 2011; Jeon et al. 2012). This is because 
accretion onto solitary BHs is more susceptible to feed- 
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back that in turn reduces the BHs accretion rate and lu- 
minosity (e.g. Johnson & Bromm 2007; Alvarez et al. 2009; 
Milosavljevic et al. 2009), whereas the near-Eddington lumi- 
nosity of a HMXB can persist through Roche lobe overflow. 
In this way the binary characteristics of Pop III stars will 
influence the nature of Pop III feedback. Also intriguing are 
observations indicating that ultra-luminous X-ray sources 
(ULXs), which may be accreting, intermediate-mass BHs, 
occur at a higher rate in less massive and lower-metallicity 
galaxies (e.g. Swartz et al. 2008). Recent modeling of binary 
system evolution over a range of metallicity also support a 
higher rate of ULX formation in lower-metallicity star clus- 
ters (Linden et al. 2010). We note that BH binaries are de- 
tectable in channels other than electromagnetic radiation. 
For instance, Kowalska et al. (2012) find that if Pop III stars 
left behind a significant BH population, then their gravita- 
tional wave (GW) emission will dominate the low-frequency 
(< 100Hz) spectrum as long as the Pop III binary fraction 
was at least 1CF 2 . 

A binary companion can also influence a star’s evolu- 
tion during and after its main sequence (MS) stage, as well as 
the type of death the star will undergo (e.g. Wellstein et al. 
2001; Petrovic et al. 2005, 2005b; Cantiello et al. 2007; 
de Mink et al. 2009; Langer 2012; Sana et al. 2012). For in- 
stance, massive stars that leave behind BH remnants may 
power collapsar gamma-ray bursts (GRBs; e.g. Woosley 
1993; Lee & Ramirez-Ruiz 2006). However, this requires suf- 
ficient angular momentum in the Pop III progenitor for 
an accretion torus to form around the remnant BH. The 
progenitor star must also lose its hydrogen envelope to 
enable the relativistic jet to penetrate through and exit 
the star (e.g. Zhang et al. 2004, but see Suwa & Ioka 2011; 
Nakauchi et al. 2012). One way to fulfill both of these condi- 
tions is for the GRB progenitor to have a close binary com- 
panion that removes the hydrogen envelope due to the heat- 
ing during a common-envelope phase (e.g. Lee et al. 2002; 
Izzard et al. 2004). 

A star which receives mass from or merges with its bi- 
nary companion potentially undergoes a process of ‘rejuve- 
nation’ in which its central hydrogen abundance increases 
(e.g. Hellings 1983; Braun & Langer 1995; Dray & Tout 
2007). This causes the star’s apparent age to be less than 
its true age, and may lead the star to appear as a blue 
straggler on the Hertzsprung- Russel diagram (Langer 2012). 
Whether rejuvenation will result from mass transfer de- 
pends on the evolutionary stage of the binary components, 
and is not expected to occur if one of the members is al- 
ready post MS. For instance, a star with an undermas- 
sive helium core may result instead (e.g. Braun & Langer 
1995). Stellar mergers may furthermore lead to the gen- 
eration of large-scale magnetic fields observed within mas- 
sive MS stars (e.g. Donati & Landstreet 2009; Ferrario et al. 
2009). In addition, while single stars must typically have 
masses > 20 Mq to go through a Wolf-Rayet (WR) phase 
(Hamann et al. 2006), and even larger masses at lower 
metallicities (Meynet & Maeder 2005), a member of a close 
binary may become a WR star at a lower mass of ~ 10 Mq 
(e.g., Vanbeveren et al. 1998; Crowther 2007). 

A star that gains mass from its binary companion 
through an accretion disk may spin up to critical rota- 
tion speed, and rapid rotation can alter a star’s evolution 
in a number of ways (see, e.g., Stacy et al. 2011, 2012; 


Maeder & Meynet 2012; Yoon et al. 2012). Rotation will al- 
low for rotationally induced mixing and possibly chemically 
homogeneous evolution (CHE). CHE, in turn, may lower the 
minimum mass at which a star will undergo a PISN to as low 
as ~ 65 Mq (e.g., Chatzopoulos & Wheeler 2012). If a ro- 
tating star instead leaves behind a BH remnant, a collapsar 
GRB may result even without a binary companion. This is 
because CHE may allow stars to bypass the red giant phase 
to become a WR star. This evolutionary path may further- 
more allow the star to retain enough angular momentum to 
become a GRB (e.g. Yoon & Langer 2005; Woosley & Heger 
2006). Rotationally induced chemical mixing will generally 
enhance a star’s luminosity, temperature, and metal pro- 
duction, though this depends sensitively on the assumed 
magnetic field evolution within the star (e.g., Ekstrom et al. 
2008; Yoon et al. 2012). Another possible consequence of 
both rotation and mass transfer is enhancement in surface 
nitrogen abundance (e.g. Ekstrom et al. 2008; Langer et al. 
2008; Brott et al. 2011; Langer 2012). 

Observations of massive stars in the Milky Way find 
that ~ 45-75% of O-stars have spectroscopic binary com- 
panions (Mason et al. 2009; Sana & Evans 2011), and calcu- 
lations of binary evolution furthermore find that 20-30% of 
O-stars will merge with their companions (Sana et al. 2012). 
Binarity therefore plays a crucial role in the evolution of 
massive stars in the Milky Way. While similar observations 
cannot be made for Pop III stars, we may still study them 
numerically. Early computational studies of Pop III star for- 
mation indicated that they would grow to be highly mas- 
sive (> 100 Mq; e.g. Abel et al. 2002; Bromm et al. 2002; 
Bromm & Loeb 2004; Yoshida et al. 2008), and that each 
minihalo would host only a single Pop III star. More recent 
work, however, has found that primordial star-forming gas 
will undergo continued fragmentation well after the initial 
star has formed. Pop III multiplicity has been found both 
in small-scale simulations that utilize idealized initial con- 
ditions (e.g. Machida et al. 2008, Clark et al 2008; 2011a), 
as well as simulations that are initialized on cosmological 
scales, (e.g. Turk et al. 2009; Stacy et al. 2010; Clark et al. 
2011b, Greif et al. 2011, 2012) and also when accounting for 
radiative feedback (e.g. Smith et al. 2011; Stacy et al. 2012). 
Recent numerical advances have thus shown that the pro- 
cesses that lead to binary formation apply to both enriched 
and metal-free gas (e.g. Tohline 2002), and a new picture has 
emerged in which Pop III stars typically form in multiples 
and within disks. 

In order to determine the role of multiplicity in Pop III 
evolution and feedback, in this study we aim to better con- 
strain the nature of Pop III binaries. We do this through 
simulating a large cosmological box which contains approx- 
imately ten minihaloes in which primordial stars form by 
2 ~ 20. We examine the rate of binary formation and pro- 
tostellar mergers, the typical characteristics of the binaries, 
and the overall mass function at this stage of evolution of 
the Pop III systems. Stars ejected from their Pop III cluster 
may stop accreting and remain low-mass and long-lived to 
the present day. Such stars may be observed in dwarf galax- 
ies or the Milky Way halo, so we furthermore discuss the fre- 
quency of low-mass sink ejections in our suite of minihaloes. 
In Section 2 we describe our numerical methods, while we 
present our results in Section 3. We address our caveats in 
Section 4, and conclude in Section 5. 
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Figure 1. Temperature versus number density of the gas just prior to the formation of the first sink particle in each minihalo. The 
similarity in the early thermal evolution is striking. At n > 10 s cm -3 , three-body reactions rapidly convert the gas into fully molecular 
form. The corresponding boost in H 2 cooling balances adiabatic heating, resulting in nearly isothermal evolution. H 2 line opacity will 
become important only at even higher densities, eventually causing deviations from the near-isothermality seen here. 


2 NUMERICAL METHODOLOGY 
2.1 Initial Setup 

We perform our investigation using gadget-2, a widely- 
tested three-dimensional N-body and SPH code (Springel 
2005). We begin with a 1.4 Mpc (comoving) box contain- 
ing 512 3 SPH gas particles and the same number of DM 
particles. The simulation is initialized at 2 = 100. Posi- 
tions and velocities are assigned to the particles in accor- 
dance with a ACDM cosmology with 12 a = 0.7, 12m = 0.3, 
Db = 0.04, as = 0.9, and h = 0.7. Each gas particle has a 
mass of m sp h = 120 Mg, while DM particles have a mass 
of mnM = 770 Mg. Once a gas particle reaches a threshold 
density of 10 3 cm" 3 , we employ ‘marker sinks’ instead of 
following these regions to increasingly higher densities. This 
allows us to efficiently evolve the entire box through a time 
period over which ten minihaloes form. Though the marker 
sinks are numerically similar to the protostellar sink parti- 
cles to be described in Section 2.4, note that the marker sinks 
are implemented at much lower densities than the protostel- 
lar sink particles of the final highest-resolution calculations. 


The location of the marker sinks coincides with the cen- 
ters of the minihaloes. In particular, we locate the first ten 
sink particles which form within the first ten minihaloes. 
Once these positions are ascertained, the cosmological box 
is reinitialized at 2 = 100 for each individual minihalo, but 
with 64 ‘child’ particles added around a 100-140 kpc region 
where the target halo will form. Larger regions of refinement 
are used for minihaloes whose mass originates from a larger 
area of the cosmological box. Particles at progressively larger 
distances from the minihalo are given increasingly larger 
masses, such that in the refined initial conditions there is 
a total of ~ 10 7 particles. The most resolved particles are of 
mass m sp h =1.85 Mq and mn m = 12 Mg. 


2.2 Cut-Out and Refinement Technique 

The refined simulations are run until the gas reaches a den- 
sity of 5 x 10' cm -3 . All particles beyond 10 physical pc 
are then removed from the simulation. Once the gas has 
reached these densities, this central star-forming region is 
gravitationally bound, and the effects of the outer minihalo 
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Figure 2. Velocity structure of the gas just prior to the formation of the first sink particle in each minihalo. Minihaloes are shown in 
order going from top left to bottom right. Solid lines show radial velocity u ra< j, dashed lines depict rotational velocity v TO t, and dash-dot 
lines show the free-fall velocity vg. Blue dotted lines represent M lur }, and follow the scaling shown on the right-hand y-axis. Turbulence 
is generally sonic (M tur |- ) =1), though in some cases the level of turbulence can become mildly supersonic. Over these distance scales, 
I;,-,,] and v TO t both remain on the order of ~ 1/2 of vg. While the H 2 chemistry leads to very similar thermal histories for each minihalo, 
differences in mass inflow history yield a greater variety in the velocity structure of the inner gas. In particular, the relative magnitude 
of u rat j and v I0 t differs in each case. Though they remain within a factor of two of each other, u ra u dominates in some cases while v ro t 
dominates in others. 


and neighboring haloes can be safely ignored. Also note that 
the gas at the cut-out edge has a free-fall time of ~ 10' yr 
and will undergo little evolution over the next 5000 yr fol- 
lowed in the simulation. Our cut-out technique leads to the 
propagation of a rarefaction wave starting from the cut-out 
edge due to the vacuum boundary condition. However, the 
wave will only travel a distance of c 3 t, where c s is the gas 
sound-speed (~ 2 kms -1 ), and the time t is 5000 yr. This 
corresponds to an insignificant distance of ~ 10 -2 pc (2000 
AU) from the cut-out edge, over two orders of magnitude 
smaller than the 10 pc box size. 

In addition, each remaining SPH particle is replaced 
with 256 child particles, each of which is placed randomly 
within the smoothing kernel of the parent particle. The mass 
of the parent particle is then evenly divided amongst the 
child particles. Each of these particles inherits the same 
chemical abundances, velocity, and entropy as the parent 
particle (see, e.g., Bromm & Loeb 2003; Clark et al. 2011b). 


This ensures conservation of mass, internal energy, and lin- 
ear momentum. After this process, each SPH particle in the 
new cut-out simulation has a mass of m ap h = 7.2 x 10 -3 Mq, 
and the new resolution mass is M les ~ 0.4 Mq, the approx- 
imate mass contained in one SPH kernel (Bate et al. 1995). 
The cut-out and refinement is performed for the central star- 
forming region within each of the ten identified minihaloes. 
We then evolve the star-forming disks in these ten different 
minihaloes over 5000 yr of protostellar accretion. 


2.3 Chemistry, Heating, and Cooling 

We utilize the same chemistry and thermal network as 
described in detail by Greif et al. (2009) and used in 
Stacy et al. (2012). Briefly, the code follows the abundance 
evolution of H, H + , H _ , H 2 , H^", He, He + , He ++ , and e~, 
as well as the three deuterium species D, D + , and HD. All 
relevant cooling mechanisms, including H 2 collisions with H 
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Figure 3. Angular momentum profile of gas within each mini- 
halo, measured just prior to the formation of the first sink. The 
thick lines of various colors and linestyles represent eight different 
minihaloes from our suite. The thin yellow lines denote angular 
momentum profiles found in separate cosmological simulations. 
Solid yellow line is taken from Stacy et al. (2010), dotted yel- 
low line from Yoshida et al. (2006), and dashed yellow line from 
Abel et al. (2002). Thick black dashed line shows an approximate 
powerlaw fit to these profiles, j oc M en c Note the near conver- 
gence of each of these angular momentum profiles, including pro- 
files obtained from varying cosmological realizations and differing 
techniques for hydrodynamics computation. 

and He as well as other H 2 molecules, are included. The 
thermal network also includes cooling through H 2 collisions 
with protons and electrons, H and He collisional excitation 
and ionization, recombination, bremsstrahlung, and inverse 
Compton scattering off the cosmic microwave background. 

Further H 2 processes are also included to properly 
model gas evolution to high densities. For instance, the 
chemistry and thermal network includes three-body H 2 for- 
mation and the concomitant H 2 formation heating, which 
become important at n > 10 8 cm~ 3 . Three-body H 2 for- 
mation rates are uncertain, however, and Turk et al. (2011) 
determined that the variation in suggested rates leads to sig- 
nificant differences in the gas collapse at high density, includ- 
ing long-term disk stability and fragmentation. For our sim- 
ulations, we choose the intermediate rate from Palla et al. 
(1983). When n > 10 9 cm -3 , cooling through H 2 ro- 
vibrational lines becomes less effective as these lines grow 
optically thick. We account for this using an escape probabil- 
ity formalism together with the Sobolev approximation (see 
Yoshida et al. 2006; Greif et al. 2011 for further details). 

2.4 Sink Particle Method 

SPH particles that reach densities of n max = 10 13 cm -3 are 
converted into sink particles. Sinks grow in mass by accret- 
ing surrounding gas particles that lie within a distance of 
r a cc, where r acc = ires — 20 AU, and 

L ies ~0.5(—') . (1) 

y P max J 


Along with the criterion that the accreted particles 
must be within a distance d < r acc from the sink, we 
also require that the particles are not rotationally sup- 
ported against infall onto the sink. This corresponds to 
jsPH < jeent, where jsph = v Iot d is the specific angular mo- 
mentum of the gas particle, j ce nt = V GM s i n kr aC c is the an- 
gular momentum required for centrifugal support, and v ro t 
and d are the rotational velocity and distance of the particle 
relative to the sink. In addition, a sink accretes all particles 
which come within d < r m i n = 4AU. These same criteria 
are used not only for gas particles but also neighboring sink 
particles, so our algorithm allows for the merging of two sink 
particles. When the sink first forms, it immediately accretes 
the majority of particles within r acc , giving it an initial mass 
of Mres — 0.4 Mq. Each time the sink grows, its position and 
velocity are set to the mass- weighted of that of the sink and 
the accreted particles. 

After a particle becomes a sink, its density, temper- 
ature, and chemical abundances are no longer updated, 
though its position and velocity continue to evolve through 
gravitational interactions. Each sink is instead held at a con- 
stant density n max and temperature of 1000 K. The sink thus 
exerts a pressure on the surrounding particles, avoiding the 
formation of an artificial pressure vacuum around its accre- 
tion radius (see Bromm et al. 2002; Martel et al. 2006). 

Once a sink is formed, nearly all gas particles near the 
sinks will be accreted before approaching the d < r m i n cri- 
terion. However, secondary sinks can come within this dis- 
tance and become subject to strong N-body interactions. 
In physical reality, however, the sub-sink region will con- 
tain not only the protostar but also surrounding gas and an 
accretion disk. In addition, the protostars themselves can 
be very distended, with photospheres reaching ~ 1 AU (e.g. 
Omukai & Palla 2003; Greif et al. 2012; Smith et al. 2012a). 
Our accretion of all secondary sinks within d < r m in there- 
fore allows us to account for sub-sink hydrodynamic affects 
which will mitigate those of pure N-body dynamics. 


2.5 Identifying Binaries 

Once a protostellar multiple system has formed, we iden- 
tify which pairs of sinks comprise a binary system. For each 
minihalo we determine which pairs of sinks are gravitation- 
ally bound by comparing their specific potential and kinetic 
energies, e p and £k- A pair of sinks is considered a binary 
system if t < 0, where e is the total orbital energy, 

e = tp + £k, (2) 

-G(Mi+M 2 ) 

£ p — r > W 

and 

£ k = in 2 . (4) 

Mi and M 2 are the masses of the two sinks in question, r is 
the distance between the sinks, and v is the relative velocity 
between the sinks. Within each minihalo we find an average 
of between two and three binaries. 

Each binary is then assigned an overall mass Mbin = 
Mi + M 2 , a center-of-mass position, and a center-of-mass 
velocity. We next determine the semi-major axis, a, for each 
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Figure 4. Density structure of the gas just prior to the formation of the first sink particle in six select minihaloes. Each star-forming 
region has a similar structure which roughly follows that of an isothermal sphere. However, the gas within differing minihaloes does vary 
somewhat in the degree of flattening. 


binary. This is derived by equating the above definition of e 
with the following energy equation: 


-GM hh 
2 a 

This yields 


( 5 ) 


—G (Mi + M 2 ) 

2 (fk + £p) 


(6) 


We furthermore calculate the binary’s mass ratio q, 
where 


Q = 


Ah 

Mi' 


( 7 ) 


In this case M 2 corresponds specifically to the less massive 
of the binary members, and Mi corresponds to the more 
massive, such that the value of q is always between zero and 
one. 


3 RESULTS 

3.1 Initial Minihalo Collapse 

The initial collapse of gas to protostellar density exhibits 
substantial similarity within each minihalo. Figure 1 shows 
the thermal state of the star-forming gas within each mini- 
halo just prior to the formation of the first high-resolution 
sink, while Figure 2 shows their velocity profiles. At high 
densities of n > 10 s cm -3 , H 2 three-body formation becomes 
rapid. This leads to nearly isothermal evolution as increased 
H 2 cooling counteracts adiabatic heating while the gas con- 
denses to n > 10 12 cm -3 . The minihaloes show a striking 


similarity in thermal evolution in this density range. In the 
evolution of gas to the densities shown in Figure 1 (n > 
10 6 cm -3 ), we also find that the HD cooling is generally 
unimportant compared to H 2 cooling, 

Figure 3 shows the angular momentum profile of the 
gas in the maximally resolved simulations just before the 
first sink has formed, measured with respect to the center- 
of-mass of all gas included in the 10 pc simulations. These 
profiles all have roughly the same magnitude and powerlaw 
slope, such that j(r) oc M en c , where M enc is the enclosed 
gas mass within a given radius r from the center-of-mass. 
This convergence is found not only among minihaloes in our 
simulation, but also including those from differing cosmolog- 
ical realizations, performed with varying methods of com- 
putational hydrodynamics (Abel et al. 2002; Yoshida et al. 
2006). These angular momentum distributions may addi- 
tionally lead to high spin rates of individual protostars (e.g. 
Stacy et al. 2011, 2012). 

The density structure of the central gas can be seen 
in Figure 4 for a representative number of the minihaloes. 
The similarity between minihaloes is again apparent, though 
some of the star-forming gas exhibits greater elongation and 
a more disk-like structure. Similar variation is seen when we 
compare the velocity and turbulent structure in Fig. 2. We 
measure M tur b over a range of radial bins according to the 
following: 

Mt ur bC 2 s = ^2 Jf ( Vi ~ Vlot ~ V ra.d) 2 , (8) 

i 

where c 3 is the sound speed of the radial bin, rm is the mass 
of a gas particle contributing to the bin, and M is the total 
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Figure 5. History of the mass growth of the stellar systems in each minihalo. Solid lines represent the total sink mass. Dashed lines 
show the growth of the most massive sink. Dotted lines depict the growth of the second-most massive sink. Average accretion rates vary 
by over a factor of three between minihaloes. 


mass of the bin. Between 20 and 10,000 AU, both the rota- 
tional velocity v ro t and radial velocity « ra d remain at a few 
km s~\ while the turbulence remains approximately sonic, 
never reaching more than twice the thermal velocity of the 
gas. In each case, both radial and infall velocities are on the 
order of ~ 1/2 of the free-fall velocity vg. However, the rela- 
tive magnitude of u ra d and v ro t varies, with u ra d dominating 
in some cases and v ro t dominating in others. 

After the central gas in each minihalo reaches the den- 
sity threshold for sink formation, in all cases a stellar mul- 
tiple system quickly forms. This fragmentation is in agree- 
ment with the parameter study of Machida et al. (2008). 
They modeled their initial conditions as modified Bonnor- 
Ebert spheres, which well-describes collapsing and approxi- 
mately isothermal gas like that within star-forming regions 
(e.g. Ebert 1955; Bonnor 1956; Stacy et al. 2010; Greif et al. 
2012), and found the following fragmentation condition: 

> Hcrit = 4 x 10~ 17 ( n ° (9) 

V lO^cm - ^ / 

where Ho and no are the angular velocity and central num- 
ber density of the initial core of the Bonnor-Ebert sphere. 
We determine Ho of the central core when no is between 


10 3 cm~ 3 and 10 4 cm -3 , corresponding to H cr i t values be- 
tween > 4 x 10- 17 s _1 and ~ 2 x 10 -16 s _1 . We take a 
mass- weighted average of the angular velocity of all gas par- 
ticles with n > 10 3 cm -3 and find in each case that Ho 
ranges between two to three orders of magnitude greater 
than Hcrit- The subsequent flattening and fragmentation of 
the dense gas is thus consistent with the criterion presented 
by Machida et al. (2008). 


3.2 Mass Function 

The mass growth history of the sinks shows substantial vari- 
ation among the ten minihaloes, as is shown in Figures 5 and 
6. At the end of the simulations, the total mass accreted by 
all sinks ranges from 30 — 120 Mg for a given minihalo (Fig. 
5). The total mass of gas within each minihalo is typically 
< 10 s Mq, implying a star-formation efficienty of ~ 10~ 3 . 
The distribution of sink masses after 5000 yr of accretion 
can be seen for the individual minihaloes in Figure 6, where 
the variation between minihaloes is even more striking. In 
particular, while some minihaloes have only one or two sinks 
with M * < 1 Mg, this number can be as large as nine. We 
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Figure 6. Distribution of sink masses in each minihalo after ~ 5000 yr of accretion. Dashed red lines depict powerlaw fits to each 
distribution, dN/dm oc M* a . We find that a ranges from zero (Minihaloes 1 and 2) to 0.57 (Minihalo 10). Dashed blue lines show an 
example distribution for a = 2, the minimum powerlaw that will yield a top-heavy mass function. Each a = 2 distribution is normalized 
to have the same total mass as that found within the sinks of the corresponding minihalo. 


find a range of powerlaw fits to the relation dN/dm oc M* a , 
where N is the number of sinks that lie within a logarith- 
mically scaled mass bin. This range extends from a = 0 for 
Minihaloes 1, 2 and 9, to a = 0.57 for Minihalo 10. Consid- 
ering the distribution across all minihaloes (Fig. 7), we find 
dN/dm oc AC 0 ' 17 . 

Unlike the observed present-day initial mass function 
(IMF), which is characterized by a Salpeter slope of a = 2.35 
towards the high mass end, our simulations yield a top-heavy 
mass function in that the majority of the protostellar mass 
resides within the most massive protostars. Top-heaviness 
applies to mass functions with a < 2, since this is the re- 
quirement for the average stellar mass to be dominated by 
the upper limit of the mass range (e.g. Bromm 2012). If we 
consider smaller values of a to be more top-heavy, then our 
fitted mass function is less top-heavy than that found after 
~ 1000 yr of sink accretion in Greif et al. (2011), who find 
a < 0 . 

Including radiative feedback would likely modify our 
values for a, though the magnitude of the modification is 
uncertain. While feedback does not prevent fragmentation 
(e.g. Smith et al. 2011, Stacy et al. 2012b), heating of the 


gas may suppress low-mass fragmentation to yield smaller 
a values closer to zero (i.e., a flatter spectrum). Greif et al. 
(2011) resolved smaller mass scales and again found smaller 
values of a. However, following such a resolved simulation 
for longer than their 1000 yr may reveal a population of very 
low-mass stars that arises at later times, instead increasing 
a. Our simulations do not resolve down to the opacity limit 
of fragmentation, the scale on which a small ~ 10~ 2 Mq 
hydrostatic core first forms (e.g. Rees 1976; Omukai & Nishi 
1998; Greif et al. 2012). We may therefore underestimate the 
number of such low-mass protostellar fragments that form. 

3.3 Binary Fraction and Merger Rate 

Every minihalo contains at least one binary, and on average 
each halo hosts between two and three binaries. We define 
the binary fraction as 


where S and B are the number of single and binary sys- 
tems, respectively. Summing over all sinks remaining in each 
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Figure 7. Total distribution of sink masses in all minihaloes. 
Red dashed line depicts a powerlaw fit to the overall distribution, 
dN/dm oc M^T 0,17 . Blue dotted lines depict the total mass dis- 
tribution from the ‘merging sink’ simulations of Greif et al. 2011 
and a powerlow fit to this distribution (a = —0.17). 

minihalo at the end of the simulations, we have S = 42 
and B = 24, yielding /b = 0.36. Calculated slightly dif- 
ferently, this means any individual star has a probability 
2 B/(S + 2B)=53% of having a binary companion. This 
is a substantial fraction, and similar to the range of mas- 
sive star binary fractions observed in various clusters in the 
Milky Way (e.g. Sana & Evans 2011; Kiminki & Kobulnicky 
2012 ). 

While a total of 217 sinks form in our suite of mini- 
haloes, 90 remain after 5000 yr. Over than half (59%) of 
the sinks are therefore lost to mergers. This is similar to 
the merger rate found in Greif et al. (2012). However, it will 
require higher-resolution simulations to determine which of 
these mergers were actually unresolved but long-lived tight 
(a < 20 AU) binaries. The results of Greif et al. (2012), 
as well as our criterion that only non-rotationally supported 
sinks are merged, support the likelihood that a large number 
of these sink mergers were indeed true protostellar mergers. 

3.4 Orbital Period and Mass Ratio Distribution 

The characteristics of the binary distribution are shown 
in Figures 8 and 9. The distribution of the binary pe- 
riods is peaked at 10 s ’ 5 days, or 870 yr. In compari- 
son with nearby solar-type stars (blue dotted line in Fig. 
8), this peak is shifted to longer periods. For instance, 
(Duquennoy & Mayor 1991) observed a P distribution that 
is peaked at 10 4 ' 8 days, or 170 years. This corresponds to 
a distribution of the semi-major axis, a, that is peaked 
around ~ 250 AU, with an average value of 760 AU. In 
Figure 8 we also show the related cumulative distribution 
of the binary periods, Fp, normalized such that Fp(P < 
oo) = 1. We display a best-fit powerlaw for this distribution, 
Fp oc log(P) 2 ' 54 , and compare to a ‘Opik distribution’ (Opik 
1925), where Fp oc log(P). Opik’s law is consistent with ob- 
servations of massive star binaries (e.g. Kouwenhoven et al. 


2007; Sana & Evans 2011; Kiminki & Kobulnicky 2012). 
However, the fit found from our simulation is significantly 
steeper. As will further be discussed in Section 4, there are 
multiple reasons for this. The spatial resolution of the simu- 
lations was limited, so the distributions cannot sample orbits 
with a < 20 AU. This contributes to the apparent lack of 
hard binaries. In addition, the binaries have only evolved 
for a few thousand years, and their orbits may evolve and 
harden over longer periods of time. 

The distribution of mass ratio q (Figure 9) goes as 
dN/dq oc g -0 ' 55 , while for the cumulative distribution we 
find a best-fit powerlaw of P q oc g 0 ' 75 . The average value of 
q is 0.35. We compare our q distribution to that found for 
solar-type stars by Duquennoy & Mayor (1991), shown as 
the blue dotted line in Fig. 9. They observed a g-distribution 
peaked at 0.23, while our distribution is peaked at g < 0.1, so 
our simulations have a comparatively high fraction of low-g 
binaries. 

Our rate of low-g binaries is also high when compared 
with more massive star clusters in the Milky Way. For in- 
stance, Kiminki & Kobulnicky (2012) find that dN/dq oc 
g 0 ' 1 , while our fit falls off much more steeply and is slightly 
outside their range of error. Sana & Evans (2011) present 
observations which yield a uniform distribution of q down 
to q = 0.2. This corresponds to the cumulative distribution 
shown as the blue line in the right panel of Fig. 9. 

If the sink particles of the simulations continued to ac- 
crete mass for longer periods of time, it is possible that 
the q distribution would shift toward a larger percentage 
of binaries with nearly equal mass. As described in, e.g., 
Bate & Bonnell (1997), if gas with sufficient angular mo- 
mentum flows onto a binary star, it will accrete preferen- 
tially onto the less massive companion. Angular momentum 
conservation prevents rotating gas from reaching the larger 
star without the aid of, e.g., gravitational torques. The gas 
instead remains at some radius beyond the more massive 
star where it is more easily captured by the less massive 
companion, resulting in a tendency for q to evolve towards 
a value of one. This is similar to the fragmentation-induced 
starvation described by Peters et al. (2010). In their simula- 
tions of star-forming molecular clouds, they found that gas 
inflowing through a stellar disk was accreted by lower-mass 
stars before it could reach the more massive and centrally- 
located star. 


3.5 Protostellar Ejection Rate 


Of the total of 90 sinks remaining at the end of our simu- 
lation, 39 have radial velocities that are greater than their 
escape velocities (Figure 10), where we define v esc as 



where M enc is the mass enclosed between the center of mass 
(COM) and the sink particle at distance r from the COM. 
The COM is determined using all gas and sink particles 
with n > 10 10 cm -3 . Roughly half (43%) of the sinks can 
thus leave the stellar disk, similar to the rate found in, e.g. 
Greif et al. (2011) in their simulations of primordial star for- 
mation and Bate et al. (2003) in their numerical studies of 
present-day, low-mass star formation. 
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Figure 8. Left: Distribution of the log of the period, P , for the binaries found within each minihalo. Note that units of P is in days. Red 
dashed line shows a Gaussian fit to the distribution. Blue dotted line shows an example fit to data for solar-type stars (Duquennoy Sz 
Mayor 1991). The distribution of the binary periods is peaked at 10 5-5 days, or 870 yr. Right: Cumulative distribution Fp of the orbital 
period of the binaries identified in all minihaloes. Red dashed line shows a powerlaw fit, Fp oc log(P) 2-54 , while the blue line shows a 
sample Opik distribution for comparison. 




Figure 9. Left: Distribution of the mass ratio q between the members of binaries located in each minihalo. Red line is the powerlaw fit 
dN/dq oc q ~ 0-55 . Blue dotted line shows an example fit to data for solar-type stars (Duquennoy Sz Mayor 1991). The average value of q is 
0.34. Right: Cumulative distribution F q of the mass ratio of all binaries in our simulation. Red dashed line is a powerlaw fit, F q oc g 0 75 . 
The blue line depicts a uniform distribution down to q = 0.2, as observed in massive star clusters in the Milky Way (e.g. Sana Sz Evans 
2011 ). 


In addition, we define a second escape velocity, that 
required to exit the minihalo: 

TW.haio = J GMhs - l ° ~ 5kms _1 , (12) 

V ^halo 

where we define Mhaio = 4 x 10 5 M 0 and rhaio = 80 pc, 
the average values for the minihaloes of our suite. Whether 
Vesc or Uesc.haio is larger varies depending on the sink and its 
location. A total of 53 sinks have u ra d > iw.haio, which is 
greater than the 39 which can escape from the stellar disk. 
Thus, the gravity of the stellar disk provides the stronger 
escape condition. In the majority of cases, for a sink to have 
sufficient velocity to escape the disk, it must already have 
sufficient velocity to escape the minihalo (left panel of Fig- 


ure 10). For a sink traveling at rw.haio, the time to cross 
the virial radius of the minihalo is ~ 2 x 10 7 yr. Note that 
within this time the minihalo is not expected to undergo sig- 
nificant mass growth itself, no more than doubling in mass 
(e.g. Stacy et al. 2011b, Wechsler et al. 2002). Thus, Uesc.haio 
will not increase by more than 1-2 km s _1 . 

Of the ejected protostars, 18 are of mass M, < 1 M 0 
and thus have a possibility of surviving to the present day 
(Clark et al 2011b). This depends, however, on whether the 
stars will accrete more mass at a later time and thus become 
shorter lived. Identifying them as Pop III stars would require 
that, as they become incorporated into the Milky Way or 
nearby dwarf galaxies, they will not accrete any enriched ma- 
terial that will remain in their atmospheres and mask them 
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Figure 10. Left: Ratio of |r ra d| to v e sc for each sink particle. Sinks with M* < 1 Mg are denoted with filled circles, while diamonds 
represent sinks with M* > 1 Mg. Horizontal dashed line shows where l^rad/^esc | = 1. Sinks above the horizontal dashed line are able to 
escape from the stellar disk. The vertical dashed line denotes v esc halo> the velocity necessary to escape the minihalo. Right: Variation of 
|f r adl with radius. Dashed line again shows v esc halo- 


as Pop II stars (e.g. Frebel et al. 2009; Johnson & Khochfar 

2011 ). 


4 CAVEATS 

The results presented in the simulation have a number of 
caveats that should be noted. First, these simulations fol- 
lowed only 5000 yr of protostellar mass accretion. The dis- 
tribution of sink masses and binary properties will continue 
to evolve beyond this point, as some of the sinks will con- 
tinue to grow in mass and orbits may harden or undergo 
disruption through close encounters with other sinks. We 
also cannot study binaries within the resolution limit of the 
simulation, 20 AU. Thus, the sink mergers discussed previ- 
ously may have instead become very tight and unresolved 
binaries. Objects formed from subsequent fragmentation on 
sub-sink scales also could not be addressed in this study. If 
we did presume that all sink mergers instead became tight 
binaries orbiting with semi- major axis of 10 AU, we would 
have a binary fraction of 73%. The average orbital radius 
of all binaries would then be 170 AU, and the average mass 
ratio would still be 0.34. 

In addition, we do not include protostellar feedback 
in our simulations. The FU-dissociating and ionizing radi- 
ation emitted from the growing protostars may substan- 
tially change the characteristics of the stellar disk and 
eventually cut off mass flow onto the disk altogether (e.g. 
Hosokawa et al. 2011; Smith et al. 2011; Stacy et al. 2012b). 
Another important physical process not included in our sim- 
ulations was magnetic fields, which may affect the angular 
momentum build-up of the protostellar disk, the subsequent 
fragmentation, and the mass flow onto the protostars (e.g. 
Tan & Blackman 2004; Silk & Langer 2006; Maki & Susa 
2007; Schleicher et al. 2010; Federrath et al. 2011; Sur et al. 
2012; Schober et al. 2012a, 2012b). In the future we plan to 


enhance the physical realism of subseqent simulations by 
including these processes and resolving smaller scales. 


5 SUMMARY AND CONCLUSIONS 

We have examined the formation of Pop III multiple sys- 
tems within ten different minihaloes extracted from a 1.4 
Mpc (comoving) cosmological box. The properties of each 
stellar group varies between minihaloes, with mass functions 
ranging from flat to a — 0.57, and average stellar accretion 
rates varying by over a factor of three between minihaloes. 

Looking at the combined properties of the ten stellar 
systems, we find that 50% of the protostars merge with 
larger protostars, while 50% of the remaining protostars are 
ejected from their stellar disk. We furthermore find a bi- 
nary fraction of ~ 36%, translating to a 53% probability 
that any given star will have a companion. Some of these 
binaries have semi-major axes that extend to 3000 AU. The 
distribution of semi-major axes peaks slightly at ~ 300 AU, 
while the mass ratio is nearly uniformly distributed between 
zero and one. Recent work by Smith et al. (2012b) implies 
that this result may be modified by the effects of heating 
and ionization from DM annihilation (DMA). If their find- 
ings hold up to further investigation, DMA may suppress 
fragmentation and lead to wider binaries of a > 1000 AU. 

Our simulation did not have sufficient resolution to re- 
solve protostellar mergers or the formation of tight binaries, 
as these processes occur on scales smaller than our resolution 
length of 20 AU. If a significant fraction of these sink merg- 
ers did represent true protostellar mergers or the formation 
of contact binaries, this has important implications for the 
Pop III stellar evolution (see, e.g. Langer 2012 and references 
therein). During a common- envelope phase, H-envelope re- 
moval by a close binary companion may allow for a star as 
low as 10 Mq to become a WR star (e.g. Vanbeveren et al. 
1998; Crowther 2007). A stellar merger may rejuvenate the 
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mass-gaining star, making it appear younger than its true 
age (e.g. Hellings 1983; Braun & Langer 1995; Dray & Tout 
2007). Alternatively, a star with an undermassive helium 
core may result (Braun & Langer 1995). Such stars may 
avoid the red supergiant phase and instead explode as 
blue supergiants. Mass accretion through Roche-lobe over- 
flow may also spin up a star, allowing for rotationally- 
induced mixing and possibly chemically homogeneous evo- 
lution, thereby altering the star’s luminosity, temperature, 
and metal production (e.g. Ekstrom et al. 2008; Yoon et al. 
2012). In particular, rotationally induced mixing may al- 
ter the nucleosynthesis within a within a low-metallicity or 
primordial massive star to yield greater CNO-element abun- 
dances, especially 14 N (e.g. Ekstrom et al. 2008; Yoon et al. 
2012). This additionally leads to a neutron excess in the 
core, enhancing the production of s-process elements (e.g. 
Pignatari et al. 2008). 

The binarity of Pop III stars also provides an important 
pathway for the formation of GRBs and high-mass X-ray bi- 
naries. If one of the binary members is in the mass range to 
collapse into a BH, then a close companion may remove its 
hydrogen envelope and allow for penetration of a GRB jet 
to the surface (e.g. Lee et al. 2002; Izzard et al. 2004). Once 
the BH remnant is left behind, the Roche lobe overflow from 
the companion enables the BH to sustain a high luminosity 
and X-ray flux. These sources could significantly contribute 
to reionization, and exert feedback effects during the as- 
sembly of the first galaxies (Mirabel et al. 2011; Jeon et al. 
2012 ). 

Pop III binarity is thus a crucial aspect to understand- 
ing Pop III evolution and feedback. A significant percent- 
age of Pop III stars exist in binaries and undergo stellar 
mergers. Greater computational power will allow for future 
simulations which include more physical processes such as 
feedback and magnetic field growth, and which will also fol- 
low the evolution of Pop III multiple systems with increased 
resolution over longer periods of time. Combined with ob- 
servational constraints, this will yield an increasingly refined 
understanding of Pop III stars and their role in the early 
Universe. 
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